Load the necessary libraries
library(car) # for regression diagnostics
library(broom) # for tidy output
library(broom.mixed)
library(ggfortify) # for model diagnostics
library(sjPlot) # for outputs
library(knitr) # for kable
library(effects) # for partial effects plots
library(ggeffects) # for effects plots in ggplot
library(emmeans) # for estimating marginal means
library(MASS) # for glm.nb
library(MuMIn) # for AICc
library(tidyverse) # for data wrangling
library(DHARMa) # for residuals and diagnostics
library(nlme) # for lme
library(lme4) # for glmer
library(glmmTMB) # for glmmTMB
library(performance) # for diagnostic plots
library(see) # for diagnostic plots
Crab_shrimp_coral
To investigate synergistic coral defence by mutualist crustaceans, Mckeon et al. (2012) conducted an aquaria experiment in which colonies of a coral species were placed in a tank along with a predatory sea star and one of four symbiont combinations:
The experiments were conducted in a large octagonal flow-through seawater tank that was partitioned into eight sections, which thereby permitted two of each of the four symbiont combinations to be observed concurrently. The tank was left overnight and in the morning, the presence of feeding scars on each coral colony was scored as evidence of predation. The experiments were repeated ten times, each time with fresh coral colonies, sea stars and symbiont.
The ten experimental times represent blocks (random effects) within which the symbiont type (fixed effect) are nested.
mckeon <- read_csv("../public/data/mckeon.csv", trim_ws = TRUE)
glimpse(mckeon)
## Rows: 80
## Columns: 3
## $ BLOCK <dbl> 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10…
## $ PREDATION <dbl> 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, …
## $ SYMBIONT <chr> "none", "none", "none", "none", "none", "none", "none", "non…
Since the response here is the presence or absence of predation (feeding scars), a binomial distribution is appropriate.
We need to make sure that the categorical variable and the random effect are defined as factors. When doing so, it might be valuable to rearrange the order of the fixed effect (SYMBIONT) such that the none group is considered the first group. This way, the other levels will all naturally be compared to this level (hence it will be treated as a reference of control group).
mckeon <- mckeon %>%
mutate(
BLOCK = factor(BLOCK),
SYMBIONT = factor(SYMBIONT, levels = c("none", "crabs", "shrimp", "both"))
)
Model formula: \[ y_i \sim{} \mathcal{N}(n, p_i)\\ ln\left(\frac{p_i}{1-p_1}\right) =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]
where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of symbionts on the probability of the colony experiencing predation. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual coral colonies.
ggplot(mckeon, aes(y = PREDATION, x = SYMBIONT)) +
geom_point(position = position_jitter(width = 0.2, height = 0)) +
facet_wrap(~BLOCK)
The data have been setup as a single factor (mixed effects) design. Alternatively, we could have coded this up as a factorial (mixed effects) such that there was a dummy variable for Crabs, and a dummy variable for Shrimps. We could then fit a model that had the main effects of Crabs and Shrimp as well as there interaction.
Since models with categorical variables essentially dummy code up the variable anyway, fitting a single factor model will result in similar parameter estimates to a factorial design with the exception of the interaction (or Crabs and Shrimp group). Nevertheless, having fit the model as a single factor model, we can then perform a specific contrast that is the same as the interaction effect from a factorial design.
Lets start with the random intercept model
mckeon.glmer1 <- glmer(PREDATION ~ SYMBIONT + (1 | BLOCK),
data = mckeon,
family = binomial(link = "logit")
)
And now attempt the random intercept/slope model. By default, glmer uses the bobyqa optimiser. This optimiser uses iterative quadratic approximation of twice differentiable functions. When functions are not twice differentiable, it may perform poorly.
We could specify the new random intercept/slope model out in full or we could simply use the update() function to make a modification (in this case, we can remove the random intercept and add the random intercept/slope).
mckeon.glmer2 <- update(mckeon.glmer1, ~ . - (1 | BLOCK) + (SYMBIONT | BLOCK))
Along with attempting to fit the model, a glmer() performs a set of model checks to help diagnose whether the optimisation process has successfully converged as well as passed other criterion.
Evidently, the above (random intercepts/slope) model failed to converge after 10000 evaluations of the likelihood function. That is, it failed to reach a consensus on the parameter estimates. Note, this is a warning not an error. It does not necessarily mean that the model fit is not correct. Nevertheless, without further exploration of what might be causing the lack of convergence, such a model should not be used as the parameter estimates (and/or their standard errors) may be unreliable.
When receiving a convergence warning, the following potential remedies are recommended:
alternative optimizers (such as Nelder-Mead, nlminbwrap or L-BFGS-B)
increase the number of iterations in an attempt to provide more scope for consensus
abandon the model
use allFit() to attempt to fit will all the available optimizers. Different optimizers perform differently under different conditions. The default optimiser (bobyqa) has been selected by the lme4 team as they consider it to be the best all-round choice and a good compromise between robustness and efficiency. However, other optimizers might work better for any given circumstance. Other common optimizers include:
If all of the optimizers yield very similar parameter estimates, then the convergence warnings can be considered as false positives and as such, ignored (unless they all suffer the same structural issue due to a poorly specified or over-fit model).
adjust the stopping tolerance for the optimiser (e.g. optCtrl=list(maxfun=2e5))
center (and possibly scale) continuous predictors
consider a more expensive Hessian calculation. Hessian computations are used for convergence checking as well as estimating the standard errors of fixed effects parameters in GLMM. By default, the method used provides very rapid estimates. However, they can be unreliable.
restart the fit with starting values taken from the current parameter estimates
Let start with the first of these options - often it is the best initial approach.
mckeon.allFit <- mckeon.glmer2 %>% allFit()
## bobyqa : [OK]
## Nelder_Mead : [OK]
## nlminbwrap : [OK]
## nmkbw : [OK]
## optimx.L-BFGS-B : [OK]
## nloptwrap.NLOPT_LN_NELDERMEAD : [OK]
## nloptwrap.NLOPT_LN_BOBYQA : [OK]
## Check which of the models are considered valid (OK)
is.OK <- sapply(mckeon.allFit, is, "merMod")
is.OK
## bobyqa Nelder_Mead
## TRUE TRUE
## nlminbwrap nmkbw
## TRUE FALSE
## optimx.L-BFGS-B nloptwrap.NLOPT_LN_NELDERMEAD
## TRUE TRUE
## nloptwrap.NLOPT_LN_BOBYQA
## TRUE
Conclusions:
We will check the convergence characteristics of only the valid models
diff_optims.OK <- mckeon.allFit[is.OK]
lapply(diff_optims.OK, function(x) x@optinfo$conv$lme4$messages)
## $bobyqa
## [1] "boundary (singular) fit: see ?isSingular"
##
## $Nelder_Mead
## [1] "boundary (singular) fit: see ?isSingular"
##
## $nlminbwrap
## [1] "boundary (singular) fit: see ?isSingular"
##
## $`optimx.L-BFGS-B`
## [1] "boundary (singular) fit: see ?isSingular"
##
## $nloptwrap.NLOPT_LN_NELDERMEAD
## [1] "unable to evaluate scaled gradient"
## [2] "Model failed to converge: degenerate Hessian with 1 negative eigenvalues"
##
## $nloptwrap.NLOPT_LN_BOBYQA
## [1] "boundary (singular) fit: see ?isSingular"
Conclusions:
Singularity refers to presence of one or more parameters whose estimated values are on the boundary of the possible parameter space (range of possible values). Typically, this is either that there are either variances that are estimated to be 0, or correlations that are estimated to be either -1 or 1.
Singular fits usually suggest over-fitting and are more likely to contain numerical problems and lead to issues concerning testing hypotheses on boundaries.
Situations that are likely to contribute to singularity include:
f|g - particularly where f is a categorical variable with a relatively large number of levels)Advice on dealing with singularity issues varies widely from simplifying the random effects structure to fitting the model in a Bayesian framework in which singularity can be avoided via appropriate priors.
For the current example, we might consider dropping the random intercept/slope.
Before excepting that we must drop the random intercept/slope, it would be worth quickly exploring the other options listed above.
Lets see if extending the number of evaluations will result in convergence.
mckeon.glmer2 <- mckeon.glmer1 %>% update(~ . - (1 | BLOCK) + (SYMBIONT | BLOCK),
control = glmerControl(
optimizer = "bobyqa",
optCtrl = list(maxfun = 2e5)
)
)
Conclusions:
In the current example, the fixed effect (SYMBIONT) is a categorical variable. Consequently, it is effectively already scaled, so the centering (and/or scaling) option isn’t a solution here.
Now lets consider the more expensive Hessian computation.
pars <- getME(mckeon.glmer2, c("theta", "fixef"))
devfun <- update(mckeon.glmer2, devFunOnly = TRUE)
library(numDeriv)
cat("hess:\n")
## hess:
print(hess <- hessian(devfun, unlist(pars)))
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.4213975722 -2.026547e-04 -1.644069e-01 1.526993e-01 -2.826724e-04
## [2,] -0.0002026547 1.473090e-05 -1.103229e-04 2.360222e-07 5.889714e-05
## [3,] -0.1644068830 -1.103229e-04 5.610199e-01 -4.732408e-01 -3.033088e-04
## [4,] 0.1526992889 2.360222e-07 -4.732408e-01 4.070909e-01 -7.901643e-05
## [5,] -0.0002826724 5.889714e-05 -3.033088e-04 -7.901643e-05 2.743109e-04
## [6,] 0.0595859600 -3.033426e-04 -1.002367e-01 1.109353e-01 -1.446770e-03
## [7,] -0.0128592645 -7.902433e-05 1.109353e-01 -8.059782e-02 -3.480331e-04
## [8,] -0.0001698081 -9.029730e-09 -1.602309e-07 -3.028307e-07 -3.996483e-09
## [9,] -0.0010049319 -2.569394e-08 -6.284820e-07 -4.679696e-07 -1.400572e-08
## [10,] -0.0009648772 -5.504486e-08 -6.568830e-07 -3.542815e-07 -1.528416e-08
## [11,] 0.0017816129 8.246683e-04 -1.147059e-02 -2.284188e-03 3.813384e-03
## [12,] 0.0005300320 -9.729221e-05 5.407133e-04 1.195055e-04 -4.498222e-04
## [13,] -0.0353598093 7.565127e-04 -3.488098e-02 -8.343249e-03 3.498610e-03
## [14,] 0.0366325791 1.654481e-04 2.286968e-02 5.939556e-03 7.655850e-04
## [,6] [,7] [,8] [,9] [,10]
## [1,] 5.958596e-02 -1.285926e-02 -1.698081e-04 -1.004932e-03 -9.648772e-04
## [2,] -3.033426e-04 -7.902433e-05 -9.029730e-09 -2.569394e-08 -5.504486e-08
## [3,] -1.002367e-01 1.109353e-01 -1.602309e-07 -6.284820e-07 -6.568830e-07
## [4,] 1.109353e-01 -8.059782e-02 -3.028307e-07 -4.679696e-07 -3.542815e-07
## [5,] -1.446770e-03 -3.480331e-04 -3.996483e-09 -1.400572e-08 -1.528416e-08
## [6,] 1.192637e-01 1.566713e-02 -1.551708e-07 -1.032991e-07 -1.750863e-07
## [7,] 1.566713e-02 5.188491e-02 9.812057e-08 -2.048719e-07 -1.002466e-07
## [8,] -1.551708e-07 9.812057e-08 5.836355e-01 -4.984751e-01 -1.088363e-03
## [9,] -1.032991e-07 -2.048719e-07 -4.984751e-01 4.260798e-01 -9.365976e-04
## [10,] -1.750863e-07 -1.002466e-07 -1.088363e-03 -9.365976e-04 4.259329e-01
## [11,] -5.303501e-02 -1.056178e-02 -9.029329e-08 -4.217722e-07 -4.492473e-07
## [12,] 2.500050e-03 5.533505e-04 -1.728246e-08 -2.747050e-08 -1.480038e-08
## [13,] -1.612714e-01 -3.857528e-02 -2.337290e-07 -6.183609e-07 -5.379129e-07
## [14,] 1.057451e-01 2.746087e-02 -7.623707e-09 -8.387474e-08 -1.079613e-07
## [,11] [,12] [,13] [,14]
## [1,] 1.781613e-03 5.300320e-04 -3.535981e-02 3.663258e-02
## [2,] 8.246683e-04 -9.729221e-05 7.565127e-04 1.654481e-04
## [3,] -1.147059e-02 5.407133e-04 -3.488098e-02 2.286968e-02
## [4,] -2.284188e-03 1.195055e-04 -8.343249e-03 5.939556e-03
## [5,] 3.813384e-03 -4.498222e-04 3.498610e-03 7.655850e-04
## [6,] -5.303501e-02 2.500050e-03 -1.612714e-01 1.057451e-01
## [7,] -1.056178e-02 5.533505e-04 -3.857528e-02 2.746087e-02
## [8,] -9.029329e-08 -1.728246e-08 -2.337290e-07 -7.623707e-09
## [9,] -4.217722e-07 -2.747050e-08 -6.183609e-07 -8.387474e-08
## [10,] -4.492473e-07 -1.480038e-08 -5.379129e-07 -1.079613e-07
## [11,] 3.960248e+00 -6.832716e-03 3.526121e-01 1.446899e-02
## [12,] -6.832716e-03 7.879640e-04 -6.203661e-03 -1.417017e-03
## [13,] 3.526121e-01 -6.203661e-03 1.164430e+00 -8.056139e-01
## [14,] 1.446899e-02 -1.417017e-03 -8.056139e-01 8.214999e-01
cat("grad:\n")
## grad:
print(grad <- grad(devfun, unlist(pars)))
## [1] 3.633557e-03 -7.390718e-09 3.913262e-08 -5.231554e-08 -4.078292e-08
## [6] 4.978622e-09 1.356960e-09 -6.296769e-08 5.468531e-08 1.552051e-08
## [11] 1.583982e-07 -2.453373e-08 -2.122501e-08 -5.751713e-08
cat("scaled gradient:\n")
## scaled gradient:
# print(scgrad <- solve(chol(hess), grad))
## compare with internal calculations:
mckeon.glmer2@optinfo$derivs
## $gradient
## [1] 3.633557e-03 -7.105427e-09 4.288125e-08 -5.375256e-08 -4.074963e-08
## [6] 4.476419e-09 1.350031e-09 -6.313172e-08 5.453415e-08 1.463718e-08
## [11] 1.515588e-07 -2.437162e-08 -2.149392e-08 -5.524470e-08
##
## $Hessian
## [,1] [,2] [,3] [,4] [,5]
## [1,] 4.206285e-01 -2.052784e-04 -1.644053e-01 1.526983e-01 -2.810955e-04
## [2,] -2.052784e-04 2.288818e-05 -1.113415e-04 7.152557e-07 5.769730e-05
## [3,] -1.644053e-01 -1.113415e-04 5.610199e-01 -4.732451e-01 -3.056526e-04
## [4,] 1.526983e-01 7.152557e-07 -4.732451e-01 4.070864e-01 -7.963181e-05
## [5,] -2.810955e-04 5.769730e-05 -3.056526e-04 -7.963181e-05 2.880096e-04
## [6,] 5.958509e-02 -3.020763e-04 -1.002369e-01 1.109383e-01 -1.446962e-03
## [7,] -1.285934e-02 -7.915497e-05 1.109328e-01 -8.060026e-02 -3.476143e-04
## [8,] 2.384186e-07 1.192093e-06 2.384186e-07 2.384186e-07 -7.152557e-07
## [9,] 1.668930e-06 -1.907349e-06 -4.768372e-07 -2.384186e-07 9.536743e-07
## [10,] 1.430511e-06 2.384186e-07 1.907349e-06 1.907349e-06 9.536743e-07
## [11,] 1.781225e-03 8.265972e-04 -1.146817e-02 -2.286196e-03 3.813505e-03
## [12,] 5.316734e-04 -9.822845e-05 5.424023e-04 1.187325e-04 -4.477501e-04
## [13,] -3.535891e-02 7.569790e-04 -3.487992e-02 -8.344889e-03 3.500223e-03
## [14,] 3.663301e-02 1.678467e-04 2.287388e-02 5.942345e-03 7.662773e-04
## [,6] [,7] [,8] [,9] [,10]
## [1,] 5.958509e-02 -1.285934e-02 2.384186e-07 1.668930e-06 1.430511e-06
## [2,] -3.020763e-04 -7.915497e-05 1.192093e-06 -1.907349e-06 2.384186e-07
## [3,] -1.002369e-01 1.109328e-01 2.384186e-07 -4.768372e-07 1.907349e-06
## [4,] 1.109383e-01 -8.060026e-02 2.384186e-07 -2.384186e-07 1.907349e-06
## [5,] -1.446962e-03 -3.476143e-04 -7.152557e-07 9.536743e-07 9.536743e-07
## [6,] 1.192703e-01 1.566839e-02 9.536743e-07 0.000000e+00 -7.152557e-07
## [7,] 1.566839e-02 5.189705e-02 1.192093e-06 -1.907349e-06 2.384186e-07
## [8,] 9.536743e-07 1.192093e-06 5.827065e-01 -4.972365e-01 -2.384186e-07
## [9,] 0.000000e+00 -1.907349e-06 -4.972365e-01 4.245358e-01 -1.668930e-06
## [10,] -7.152557e-07 2.384186e-07 -2.384186e-07 -1.668930e-06 4.245272e-01
## [11,] -5.303502e-02 -1.056194e-02 -4.768372e-07 -4.768372e-07 -4.768372e-07
## [12,] 2.501011e-03 5.519390e-04 -2.384186e-06 2.622604e-06 7.152557e-07
## [13,] -1.612737e-01 -3.857756e-02 -2.384186e-07 2.384186e-06 1.192093e-06
## [14,] 1.057370e-01 2.746010e-02 -1.907349e-06 2.145767e-06 -4.768372e-07
## [,11] [,12] [,13] [,14]
## [1,] 1.781225e-03 5.316734e-04 -3.535891e-02 3.663301e-02
## [2,] 8.265972e-04 -9.822845e-05 7.569790e-04 1.678467e-04
## [3,] -1.146817e-02 5.424023e-04 -3.487992e-02 2.287388e-02
## [4,] -2.286196e-03 1.187325e-04 -8.344889e-03 5.942345e-03
## [5,] 3.813505e-03 -4.477501e-04 3.500223e-03 7.662773e-04
## [6,] -5.303502e-02 2.501011e-03 -1.612737e-01 1.057370e-01
## [7,] -1.056194e-02 5.519390e-04 -3.857756e-02 2.746010e-02
## [8,] -4.768372e-07 -2.384186e-06 -2.384186e-07 -1.907349e-06
## [9,] -4.768372e-07 2.622604e-06 2.384186e-06 2.145767e-06
## [10,] -4.768372e-07 7.152557e-07 1.192093e-06 -4.768372e-07
## [11,] 3.960271e+00 -6.834507e-03 3.526089e-01 1.446557e-02
## [12,] -6.834507e-03 7.982254e-04 -6.203413e-03 -1.418829e-03
## [13,] 3.526089e-01 -6.203413e-03 1.164444e+00 -8.056104e-01
## [14,] 1.446557e-02 -1.418829e-03 -8.056104e-01 8.215151e-01
Conclusions:
Finally, lets consider re-fitting the model using the current parameter estimates as starting values.
# strict_tol <- glmerControl(optCtrl=list(xtol_abs=1e-8, ftol_abs=1e-8))
pars <- getME(mckeon.glmer2, c("theta", "fixef"))
mckeon.glmer3 <- update(mckeon.glmer2, start = pars, control = glmerControl(
optimizer = "bobyqa",
optCtrl = list(maxfun = 2e5)
))
Conclusions:
Lets start with the random intercept model
mckeon.glmmTMB1 <- glmmTMB(PREDATION ~ SYMBIONT + (1 | BLOCK),
data = mckeon,
family = binomial(link = "logit"),
REML = TRUE
)
And now attempt the random intercept/slope model.
mckeon.glmmTMB2 <- mckeon.glmmTMB1 %>% update(~ . - (1 | BLOCK) + (SYMBIONT | BLOCK))
## Error in (function (start, objective, gradient = NULL, hessian = NULL, : NA/NaN gradient evaluation
mckeon.glmmTMB2 <- mckeon.glmmTMB1 %>% update(~ . - (1 | BLOCK) + (SYMBIONT | BLOCK),
control = glmmTMBControl(
optimizer = optim,
optArgs = list(method = "SANN")
)
)
AICc(mckeon.glmmTMB1, mckeon.glmmTMB2)
IN addition to numerous warnings, there is an error about NA/NaN gradient evaluation. This is a somewhat generic error returned by the nlminb optimiser. It occurs when there is an error in either the objective function or gradient function calculations.
Ultimately, this might suggest a poorly specified or over-fit model. We will drop the random intercept/slope.
mckeon.glmer1 %>%
plot_model(type = "diag") %>%
plot_grid()
Conclusions:
mckeon.glmer1 %>% performance::check_model()
Conclusions:
mckeon.resid <- mckeon.glmer1 %>% simulateResiduals(plot = TRUE)
mckeon.resid %>% testZeroInflation()
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0428, p-value = 0.936
## alternative hypothesis: two.sided
Conclusions:
mckeon.glmmTMB1 %>%
plot_model(type = "diag") %>%
plot_grid()
Conclusions:
mckeon.glmmTMB1 %>% performance::check_model()
Conclusions:
mckeon.resid <- mckeon.glmmTMB1 %>% simulateResiduals(plot = TRUE)
mckeon.resid %>% testZeroInflation()
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.96105, p-value = 0.968
## alternative hypothesis: two.sided
mckeon.resid %>% testDispersion()
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 1.0036, p-value = 0.928
## alternative hypothesis: two.sided
Conclusions:
mckeon.glmer1 %>% plot_model(type = "eff", show.data = TRUE, jitter = c(0.05, 0))
## $SYMBIONT
mckeon.glmer1 %>%
allEffects() %>%
plot()
mckeon.glmer1 %>%
ggpredict() %>%
plot(add.data = TRUE, jitter = c(0.1, 0))
## $SYMBIONT
mckeon.glmer1 %>%
ggemmeans(~SYMBIONT) %>%
plot(add.data = TRUE, jitter = c(0.1, 0))
mckeon.glmmTMB1 %>% plot_model(type = "eff", show.data = TRUE, jitter = c(0.05, 0))
## $SYMBIONT
mckeon.glmmTMB1 %>%
allEffects() %>%
plot()
mckeon.glmmTMB1 %>%
ggpredict() %>%
plot(add.data = TRUE, jitter = c(0.05, 0))
## $SYMBIONT
mckeon.glmmTMB1 %>%
ggemmeans(~SYMBIONT) %>%
plot(add.data = TRUE, jitter = c(0.05, 0))
mckeon.glmer1 %>% summary()
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: PREDATION ~ SYMBIONT + (1 | BLOCK)
## Data: mckeon
##
## AIC BIC logLik deviance df.resid
## 70.7 82.6 -30.4 60.7 75
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -23.2076 -0.1730 0.0976 0.2980 0.8944
##
## Random effects:
## Groups Name Variance Std.Dev.
## BLOCK (Intercept) 11.81 3.437
## Number of obs: 80, groups: BLOCK, 10
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.096 1.811 2.814 0.00490 **
## SYMBIONTcrabs -3.842 1.465 -2.623 0.00871 **
## SYMBIONTshrimp -4.431 1.551 -2.856 0.00429 **
## SYMBIONTboth -5.599 1.724 -3.247 0.00116 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) SYMBIONTc SYMBIONTs
## SYMBIONTcrb -0.625
## SYMBIONTshr -0.652 0.736
## SYMBIONTbth -0.687 0.732 0.767
Conclusions:
mckeon.glmer1 %>% tidy(effect = "fixed", conf.int = TRUE, infer = c(TRUE, TRUE), exponentiate = TRUE)
mckeon.glmer1 %>%
tidy(effect = "fixed", conf.int = TRUE, infer = c(TRUE, TRUE), exponentiate = TRUE) %>%
kable()
| effect | term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|---|
| fixed | (Intercept) | 163.3256181 | 295.7986843 | 2.813623 | 0.0048987 | 4.6929398 | 5684.1252579 |
| fixed | SYMBIONTcrabs | 0.0214459 | 0.0314110 | -2.623280 | 0.0087088 | 0.0012151 | 0.3785015 |
| fixed | SYMBIONTshrimp | 0.0119023 | 0.0184657 | -2.856063 | 0.0042893 | 0.0005689 | 0.2490134 |
| fixed | SYMBIONTboth | 0.0037016 | 0.0063822 | -3.247315 | 0.0011650 | 0.0001261 | 0.1086478 |
Conclusions:
# warning this is only appropriate for html output
mckeon.glmer1 %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)
| PREDATION | ||||
|---|---|---|---|---|
| Predictors | Odds Ratios | std. Error | CI | p |
| (Intercept) | 163.33 | 295.80 | 4.69 – 5684.13 | 0.005 |
| SYMBIONT [crabs] | 0.02 | 0.03 | 0.00 – 0.38 | 0.009 |
| SYMBIONT [shrimp] | 0.01 | 0.02 | 0.00 – 0.25 | 0.004 |
| SYMBIONT [both] | 0.00 | 0.01 | 0.00 – 0.11 | 0.001 |
| Random Effects | ||||
| σ2 | 3.29 | |||
| τ00 BLOCK | 11.81 | |||
| ICC | 0.78 | |||
| N BLOCK | 10 | |||
| Observations | 80 | |||
| Marginal R2 / Conditional R2 | 0.228 / 0.832 | |||
| AIC | 70.706 | |||
mckeon.glmmTMB1 %>% summary()
## Family: binomial ( logit )
## Formula: PREDATION ~ SYMBIONT + (1 | BLOCK)
## Data: mckeon
##
## AIC BIC logLik deviance df.resid
## 62.9 74.8 -26.4 52.9 79
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## BLOCK (Intercept) 17.83 4.223
## Number of obs: 80, groups: BLOCK, 10
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.420 1.782 2.481 0.013109 *
## SYMBIONTcrabs -3.317 1.313 -2.526 0.011534 *
## SYMBIONTshrimp -3.956 1.416 -2.793 0.005218 **
## SYMBIONTboth -5.152 1.565 -3.291 0.000997 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Conclusions:
mckeon.glmmTMB1 %>% tidy(effect = "fixed", conf.int = TRUE, infer = c(TRUE, TRUE))
mckeon.glmmTMB1 %>% tidy(effect = "fixed", conf.int = TRUE, infer = c(TRUE, TRUE), exponentiate = TRUE)
Conclusions:
# warning this is only appropriate for html output
mckeon.glmmTMB1 %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)
| PREDATION | ||||
|---|---|---|---|---|
| Predictors | Odds Ratios | std. Error | CI | p |
| (Intercept) | 83.10 | 148.05 | 2.53 – 2730.03 | 0.013 |
| SYMBIONT [crabs] | 0.04 | 0.05 | 0.00 – 0.48 | 0.012 |
| SYMBIONT [shrimp] | 0.02 | 0.03 | 0.00 – 0.31 | 0.005 |
| SYMBIONT [both] | 0.01 | 0.01 | 0.00 – 0.12 | 0.001 |
| Random Effects | ||||
| σ2 | 3.29 | |||
| τ00 BLOCK | 17.83 | |||
| ICC | 0.84 | |||
| N BLOCK | 10 | |||
| Observations | 80 | |||
| Marginal R2 / Conditional R2 | 0.149 / 0.867 | |||
| AIC | 62.861 | |||
In addition to comparing each of the symbiont types against the control group of no symbionts, it might be interesting to investigate whether there are any differences between the predation protection provided by crabs and shrimp, as well as whether having both crabs and shrimp symbionts is different to only a single symbiont type.
These contrasts can be explored via specific contrasts.
| SYMBIONT | Crab vs Shrimp | One vs Both | None vs Symbiont |
|---|---|---|---|
| none | 0 | 0 | 1 |
| crab | 1 | 1/2 | -1/3 |
| shrimp | -1 | 1/2 | -1/3 |
| both | 0 | -1 | -1/3 |
cmat <- cbind(
crab_vs_shrimp = c(0, 1, -1, 0),
one_vs_both = c(0, 1 / 2, 1 / 2, -1),
symbiont = c(1, -1 / 3, -1 / 3, -1 / 3)
)
round(crossprod(cmat), 1)
## crab_vs_shrimp one_vs_both symbiont
## crab_vs_shrimp 2 0.0 0.0
## one_vs_both 0 1.5 0.0
## symbiont 0 0.0 1.3
## all contrasts orthogonal
mckeon.glmer1 %>%
emmeans(~SYMBIONT, type = "response") %>%
contrast(method = list(SYMBIONT = cmat)) %>%
summary(infer = TRUE)
## Or on an absolute difference scale
mckeon.glmer1 %>%
emmeans(~SYMBIONT, link = "link") %>%
regrid() %>%
contrast(method = list(SYMBIONT = cmat)) %>%
summary(infer = TRUE)
Conclusions:
mckeon.glmer1 %>% r.squaredGLMM()
## R2m R2c
## theoretical 0.2282003 0.8318487
## delta 0.2166278 0.7896639
## The delta mehod can be used with for all distributions and link
## functions, jwhile lognormal approximation and trigamma function are
## limited to distributions with logarithmic link. Trigamma-estimate
## is recommended whenever available. Additionally, for binomial
## distributions, theoretical variances exist specific for each link
## function distribution.
mckeon.glmer1 %>% performance::r2_nakagawa()
## # R2 for Mixed Models
##
## Conditional R2: 0.832
## Marginal R2: 0.228
cmat <- cbind(
crab_vs_shrimp = c(0, 1, -1, 0),
one_vs_both = c(0, 1 / 2, 1 / 2, -1),
symbiont = c(1, -1 / 3, -1 / 3, -1 / 3)
)
round(crossprod(cmat), 1)
## crab_vs_shrimp one_vs_both symbiont
## crab_vs_shrimp 2 0.0 0.0
## one_vs_both 0 1.5 0.0
## symbiont 0 0.0 1.3
## all contrasts orthogonal
## On the link scale
mckeon.glmmTMB1 %>%
emmeans(~SYMBIONT, type = "link") %>%
contrast(method = list(SYMBIONT = cmat)) %>%
summary(infer = TRUE)
## On the ratio (fold) scale
mckeon.glmmTMB1 %>%
emmeans(~SYMBIONT, type = "response") %>%
contrast(method = list(SYMBIONT = cmat)) %>%
summary(infer = TRUE)
## On an absolute difference scale
mckeon.glmmTMB1 %>%
emmeans(~SYMBIONT, type = "link") %>%
regrid() %>%
contrast(method = list(SYMBIONT = cmat)) %>%
summary(infer = TRUE)
Conclusions:
mckeon.glmmTMB1 %>% r.squaredGLMM()
## R2m R2c
## theoretical 0.1489321 0.8674549
## delta 0.1437647 0.8373578
## The delta method can be used with for all distributions and link
## functions, while lognormal approximation and trigamma function are
## limited to distributions with logarithmic link. Trigamma-estimate
## is recommended whenever available. Additionally, for binomial
## distributions, theoretical variances exist specific for each link
## function distribution.
mckeon.glmmTMB1 %>% performance::r2_nakagawa()
## # R2 for Mixed Models
##
## Conditional R2: 0.867
## Marginal R2: 0.149
mckeon.glmer1 %>%
emmeans(~SYMBIONT, type = "response") %>%
as.data.frame() %>%
ggplot(aes(y = prob, x = SYMBIONT)) +
geom_pointrange(aes(ymin = asymp.LCL, ymax = asymp.UCL))
mckeon.glmmTMB1 %>%
emmeans(~SYMBIONT, type = "response") %>%
as.data.frame() %>%
ggplot(aes(y = prob, x = SYMBIONT)) +
geom_pointrange(aes(ymin = lower.CL, ymax = upper.CL))